"""
Comprehensive tensor field implementation for PDEJulia.jl

Based on py-pde tensorial field implementation with Julia optimizations.
Implements rank-2 tensor fields on grids with full tensor operations.
"""

using LinearAlgebra
using StaticArrays
using Statistics
using ..Grids

"""
    Tensor2Field{T<:Real, N, G<:AbstractGrid{N}} <: AbstractDataField

A rank-2 tensor field defined on an N-dimensional grid.

The data is stored as an (N+2)-dimensional array where the last two dimensions
correspond to tensor indices. For a 2D grid, the data shape is (nx, ny, 2, 2),
where data[:, :, i, j] is the (i,j) component of the tensor at each grid point.

Warning:
    Components of the tensor field are given in the local basis. While the local
    basis is identical to the global basis in Cartesian coordinates, the local
    basis depends on position in curvilinear coordinate systems.

# Fields
- `_grid::G`: The underlying spatial grid
- `_data::Array{T,N+2}`: The field data array with shape (grid_dims..., N, N)

# Examples
```julia
# Create from function  
grid = UnitGrid([32, 32])
field = Tensor2Field(pos -> [pos[1] 0; 0 pos[2]], grid)

# Create identity tensor field
field = Tensor2Field(grid) do pos
    Matrix{Float64}(I, 2, 2)
end

# Create zero tensor field
field = Tensor2Field(Float64, grid)
```
"""
struct Tensor2Field{T<:Real, N, G<:AbstractGrid{N}} <: AbstractDataField
    _grid::G
    _data::Array{T,N+2}  # Shape: (grid_dims..., N, N)
    
    function Tensor2Field{T,N,G}(grid::G, data::Array{T,N+2}) where {T<:Real, N, G<:AbstractGrid{N}}
        expected_size = (size(grid)..., N, N)
        if size(data) != expected_size
            throw(DimensionMismatch(
                "Data size $(size(data)) does not match expected $(expected_size)"
            ))
        end
        new{T,N,G}(grid, data)
    end
end

# Interface implementation  
grid(field::Tensor2Field) = field._grid
data(field::Tensor2Field) = field._data

# Constructors

"""
    Tensor2Field(data::Array{T,N+2}, grid::AbstractGrid{N}) where {T<:Real, N}

Create tensor field from data array and grid.
"""
function Tensor2Field(data::Array{T,N+2}, grid::G) where {T<:Real, N, G<:AbstractGrid{N}}
    Tensor2Field{T,N,G}(grid, copy(data))
end

"""
    Tensor2Field(func::Function, grid::AbstractGrid{N}, ::Type{T}=Float64) where {T<:Real, N}

Create tensor field by evaluating function at grid points.
The function should take a coordinate vector and return an N×N matrix.
"""
function Tensor2Field(func::Function, grid::G, ::Type{T}=Float64) where {T<:Real, N, G<:AbstractGrid{N}}
    data = Array{T,N+2}(undef, size(grid)..., N, N)
    
    for idx in CartesianIndices(size(grid))
        coord = coordinate(grid, Tuple(idx))
        tensor_val = func(coord)
        
        # Handle both Matrix and static array inputs
        if tensor_val isa AbstractMatrix
            for i in 1:N, j in 1:N
                data[idx, i, j] = T(tensor_val[i, j])
            end
        else
            throw(ArgumentError("Function must return an N×N matrix"))
        end
    end
    
    Tensor2Field{T,N,G}(grid, data)
end

"""
    Tensor2Field(::Type{T}, grid::AbstractGrid{N}) where {T<:Real, N}

Create zero-initialized tensor field of given type.
"""
function Tensor2Field(::Type{T}, grid::G) where {T<:Real, N, G<:AbstractGrid{N}}
    data = zeros(T, size(grid)..., N, N)
    Tensor2Field{T,N,G}(grid, data)
end

"""
    Tensor2Field(grid::AbstractGrid{N}) where N

Create zero-initialized tensor field with Float64 type.
"""
Tensor2Field(grid::AbstractGrid) = Tensor2Field(Float64, grid)

"""
    identity_tensor(grid::AbstractGrid{N}, ::Type{T}=Float64) where {T<:Real, N}

Create identity tensor field (I at each grid point).
"""
function identity_tensor(grid::G, ::Type{T}=Float64) where {T<:Real, N, G<:AbstractGrid{N}}
    data = zeros(T, size(grid)..., N, N)
    
    # Set diagonal elements to 1
    for idx in CartesianIndices(size(grid))
        for i in 1:N
            data[idx, i, i] = one(T)
        end
    end
    
    Tensor2Field{T,N,G}(grid, data)
end

# Base interface
function Base.similar(field::Tensor2Field{T,N,G}) where {T,N,G}
    Tensor2Field{T,N,G}(field._grid, similar(field._data))
end

function Base.copy(field::Tensor2Field{T,N,G}) where {T,N,G}
    Tensor2Field{T,N,G}(field._grid, copy(field._data))
end

Base.size(field::Tensor2Field) = size(field._grid)
Base.length(field::Tensor2Field) = length(field._grid) * ndim(field._grid)^2
Base.eltype(field::Tensor2Field{T}) where T = T
Base.ndims(field::Tensor2Field{T,N}) where {T,N} = N

# Component access

"""
    getindex(field::Tensor2Field, i::Int, j::Int) -> ScalarField
    field[i,j]

Get the (i,j) component as a scalar field.
"""
function Base.getindex(field::Tensor2Field{T,N,G}, i::Int, j::Int) where {T,N,G}
    if !(1 ≤ i ≤ N) || !(1 ≤ j ≤ N)
        throw(BoundsError(field, (i, j)))
    end
    
    component_data = data(field)[ntuple(d -> d == N+1 ? i : (d == N+2 ? j : Colon()), N+2)...]
    ScalarField{T}(grid(field), component_data)
end

"""
    setindex!(field::Tensor2Field, value, i::Int, j::Int)
    field[i,j] = value

Set the (i,j) component from a scalar field or numeric value.
"""
function Base.setindex!(field::Tensor2Field{T,N}, value::ScalarField{T}, i::Int, j::Int) where {T,N}
    if !(1 ≤ i ≤ N) || !(1 ≤ j ≤ N)
        throw(BoundsError(field, (i, j)))
    end
    
    compatible_fields(field[1,1], value) || throw(ArgumentError("Incompatible grids"))
    
    data(field)[ntuple(d -> d == N+1 ? i : (d == N+2 ? j : Colon()), N+2)...] .= data(value)
end

function Base.setindex!(field::Tensor2Field{T,N}, value::Number, i::Int, j::Int) where {T,N}
    if !(1 ≤ i ≤ N) || !(1 ≤ j ≤ N)
        throw(BoundsError(field, (i, j)))
    end
    
    data(field)[ntuple(d -> d == N+1 ? i : (d == N+2 ? j : Colon()), N+2)...] .= T(value)
end

# Matrix extraction at specific points
"""
    get_matrix(field::Tensor2Field{T,N}, idx...) -> Matrix{T}

Extract the N×N tensor matrix at the given grid indices.
"""
function get_matrix(field::Tensor2Field{T,N}, idx...) where {T,N}
    matrix = Matrix{T}(undef, N, N)
    for i in 1:N, j in 1:N
        matrix[i,j] = data(field)[idx..., i, j]
    end
    matrix
end

"""
    set_matrix!(field::Tensor2Field{T,N}, matrix::Matrix{T}, idx...)

Set the N×N tensor matrix at the given grid indices.
"""
function set_matrix!(field::Tensor2Field{T,N}, matrix::Matrix{T}, idx...) where {T,N}
    if size(matrix) != (N, N)
        throw(DimensionMismatch("Matrix must be $(N)×$(N), got $(size(matrix))"))
    end
    
    for i in 1:N, j in 1:N
        data(field)[idx..., i, j] = matrix[i,j]
    end
end# Extended Tensor2Field functionality

# Tensor operations

"""
    trace(field::Tensor2Field) -> ScalarField

Compute the trace of the tensor field at each grid point.
"""
function LinearAlgebra.trace(field::Tensor2Field{T,N}) where {T,N}
    result_data = zeros(T, size(grid(field)))
    
    for i in 1:N
        result_data .+= data(field)[ntuple(d -> d == N+1 ? i : (d == N+2 ? i : Colon()), N+2)...]
    end
    
    ScalarField{T}(grid(field), result_data)
end

"""
    determinant(field::Tensor2Field) -> ScalarField

Compute the determinant of the tensor field at each grid point.
"""
function LinearAlgebra.det(field::Tensor2Field{T,N}) where {T,N}
    result_data = Array{T}(undef, size(grid(field)))
    
    if N == 1
        # 1×1 case: determinant is just the single element
        result_data .= data(field)[:, 1, 1]
    elseif N == 2
        # 2×2 case: ad - bc
        a = data(field)[ntuple(d -> d == N+1 ? 1 : (d == N+2 ? 1 : Colon()), N+2)...]
        b = data(field)[ntuple(d -> d == N+1 ? 1 : (d == N+2 ? 2 : Colon()), N+2)...]
        c = data(field)[ntuple(d -> d == N+1 ? 2 : (d == N+2 ? 1 : Colon()), N+2)...]
        d = data(field)[ntuple(d -> d == N+1 ? 2 : (d == N+2 ? 2 : Colon()), N+2)...]
        result_data .= a .* d .- b .* c
    elseif N == 3
        # 3×3 case: full determinant formula
        for idx in CartesianIndices(size(grid(field)))
            matrix = get_matrix(field, Tuple(idx)...)
            result_data[idx] = LinearAlgebra.det(matrix)
        end
    else
        # General case: use LinearAlgebra.det for each matrix
        for idx in CartesianIndices(size(grid(field)))
            matrix = get_matrix(field, Tuple(idx)...)
            result_data[idx] = LinearAlgebra.det(matrix)
        end
    end
    
    ScalarField{T}(grid(field), result_data)
end

"""
    transpose(field::Tensor2Field) -> Tensor2Field

Compute the transpose of the tensor field.
"""
function LinearAlgebra.transpose(field::Tensor2Field{T,N,G}) where {T,N,G}
    result_data = Array{T,N+2}(undef, size(data(field)))
    
    for i in 1:N, j in 1:N
        # Transpose: result[i,j] = field[j,i]
        result_data[ntuple(d -> d == N+1 ? i : (d == N+2 ? j : Colon()), N+2)...] =
            data(field)[ntuple(d -> d == N+1 ? j : (d == N+2 ? i : Colon()), N+2)...]
    end
    
    Tensor2Field{T,N,G}(grid(field), result_data)
end

"""
    is_symmetric(field::Tensor2Field; rtol=1e-5, atol=1e-8) -> Bool

Check if the tensor field is symmetric at all grid points.
"""
function is_symmetric(field::Tensor2Field{T,N}; rtol=1e-5, atol=1e-8) where {T,N}
    field_T = transpose(field)
    return isapprox(data(field), data(field_T); rtol=rtol, atol=atol)
end

"""
    symmetrize(field::Tensor2Field; make_traceless=false) -> Tensor2Field

Symmetrize the tensor field: (T + T^T)/2. Optionally make it traceless.
"""
function symmetrize(field::Tensor2Field{T,N,G}; make_traceless=false) where {T,N,G}
    result = (field + transpose(field)) / 2
    
    if make_traceless
        tr = trace(result)
        trace_part = tr / N
        
        # Subtract trace part from diagonal
        for i in 1:N
            result[i,i] = result[i,i] - trace_part
        end
    end
    
    result
end

"""
    eigenvalues(field::Tensor2Field) -> Vector{ScalarField}

Compute eigenvalues of the tensor field at each grid point.
Returns a vector of scalar fields, one for each eigenvalue.
"""
function LinearAlgebra.eigenvals(field::Tensor2Field{T,N}) where {T,N}
    eigenvalue_fields = [ScalarField{T}(grid(field)) for _ in 1:N]
    
    for idx in CartesianIndices(size(grid(field)))
        matrix = get_matrix(field, Tuple(idx)...)
        eigvals = LinearAlgebra.eigenvals(matrix)
        
        # Sort eigenvalues for consistency
        sort!(eigvals, rev=true)  # Descending order
        
        for (i, λ) in enumerate(eigvals)
            data(eigenvalue_fields[i])[idx] = real(λ)  # Take real part
        end
    end
    
    eigenvalue_fields
end

"""
    eigenvectors(field::Tensor2Field) -> Tuple{Vector{ScalarField}, Vector{VectorField}}

Compute eigenvalues and eigenvectors of the tensor field.
Returns (eigenvalues, eigenvectors) where eigenvectors[i] corresponds to eigenvalues[i].
"""
function LinearAlgebra.eigen(field::Tensor2Field{T,N,G}) where {T,N,G}
    eigenvalue_fields = [ScalarField{T}(grid(field)) for _ in 1:N]
    eigenvector_fields = [VectorField{T,N,G}(grid(field)) for _ in 1:N]
    
    for idx in CartesianIndices(size(grid(field)))
        matrix = get_matrix(field, Tuple(idx)...)
        
        # Handle potential complex eigenvalues for non-symmetric matrices
        if is_symmetric_matrix(matrix)
            # Use specialized symmetric eigen for better accuracy
            eigvals, eigvecs = LinearAlgebra.eigen(Symmetric(matrix))
        else
            eigvals, eigvecs = LinearAlgebra.eigen(matrix)
        end
        
        # Sort by eigenvalue magnitude (descending)
        perm = sortperm(real.(eigvals), rev=true)
        eigvals = eigvals[perm]
        eigvecs = eigvecs[:, perm]
        
        for i in 1:N
            data(eigenvalue_fields[i])[idx] = real(eigvals[i])
            
            # Store eigenvector components
            for j in 1:N
                eigenvector_data = data(eigenvector_fields[i])
                eigenvector_data[Tuple(idx)..., j] = real(eigvecs[j, i])
            end
        end
    end
    
    eigenvalue_fields, eigenvector_fields
end

# Helper function for checking matrix symmetry
function is_symmetric_matrix(matrix::AbstractMatrix; rtol=1e-12)
    return isapprox(matrix, transpose(matrix); rtol=rtol)
end

# Tensor invariants

"""
    first_invariant(field::Tensor2Field) -> ScalarField

Compute the first invariant I₁ = tr(T).
"""
first_invariant(field::Tensor2Field) = trace(field)

"""
    second_invariant(field::Tensor2Field) -> ScalarField

Compute the second invariant I₂ = ½[tr(T)² - tr(T²)].
"""
function second_invariant(field::Tensor2Field{T,N}) where {T,N}
    result_data = zeros(T, size(grid(field)))
    
    for idx in CartesianIndices(size(grid(field)))
        matrix = get_matrix(field, Tuple(idx)...)
        tr_T = LinearAlgebra.tr(matrix)
        tr_T2 = LinearAlgebra.tr(matrix * matrix)
        result_data[idx] = 0.5 * (tr_T^2 - tr_T2)
    end
    
    ScalarField{T}(grid(field), result_data)
end

"""
    third_invariant(field::Tensor2Field) -> ScalarField

Compute the third invariant I₃ = det(T).
"""
third_invariant(field::Tensor2Field) = det(field)

# Tensor products and contractions

"""
    dot(tensor::Tensor2Field, vector::VectorField) -> VectorField
    tensor ⋅ vector

Compute tensor-vector product: (T⋅v)ᵢ = Σⱼ Tᵢⱼvⱼ.
"""
function LinearAlgebra.dot(tensor::Tensor2Field{T1,N,G}, vector::VectorField{T2,N,G}) where {T1,T2,N,G}
    compatible_fields(tensor[1,1], vector[1]) || throw(ArgumentError("Incompatible grids"))
    
    T = promote_type(T1, T2)
    result_data = Array{T,N+1}(undef, size(grid(tensor))..., N)
    
    for i in 1:N
        component_data = zeros(T, size(grid(tensor)))
        for j in 1:N
            tensor_comp = data(tensor)[ntuple(d -> d == N+1 ? i : (d == N+2 ? j : Colon()), N+2)...]
            vector_comp = data(vector)[ntuple(d -> d == N+1 ? j : Colon(), N+1)...]
            component_data .+= tensor_comp .* vector_comp
        end
        result_data[ntuple(d -> d == N+1 ? i : Colon(), N+1)...] = component_data
    end
    
    VectorField{T,N,G}(grid(tensor), result_data)
end

"""
    dot(vector::VectorField, tensor::Tensor2Field) -> VectorField
    vector ⋅ tensor

Compute vector-tensor product: (v⋅T)ⱼ = Σᵢ vᵢTᵢⱼ.
"""
function LinearAlgebra.dot(vector::VectorField{T1,N,G}, tensor::Tensor2Field{T2,N,G}) where {T1,T2,N,G}
    compatible_fields(vector[1], tensor[1,1]) || throw(ArgumentError("Incompatible grids"))
    
    T = promote_type(T1, T2)
    result_data = Array{T,N+1}(undef, size(grid(tensor))..., N)
    
    for j in 1:N
        component_data = zeros(T, size(grid(tensor)))
        for i in 1:N
            vector_comp = data(vector)[ntuple(d -> d == N+1 ? i : Colon(), N+1)...]
            tensor_comp = data(tensor)[ntuple(d -> d == N+1 ? i : (d == N+2 ? j : Colon()), N+2)...]
            component_data .+= vector_comp .* tensor_comp
        end
        result_data[ntuple(d -> d == N+1 ? j : Colon(), N+1)...] = component_data
    end
    
    VectorField{T,N,G}(grid(tensor), result_data)
end

"""
    double_dot(t1::Tensor2Field, t2::Tensor2Field) -> ScalarField
    t1 : t2

Compute double contraction (Frobenius inner product): Σᵢⱼ T1ᵢⱼT2ᵢⱼ.
"""
function double_dot(t1::Tensor2Field{T1,N,G}, t2::Tensor2Field{T2,N,G}) where {T1,T2,N,G}
    compatible_fields(t1[1,1], t2[1,1]) || throw(ArgumentError("Incompatible grids"))
    
    T = promote_type(T1, T2)
    result_data = zeros(T, size(grid(t1)))
    
    for i in 1:N, j in 1:N
        t1_comp = data(t1)[ntuple(d -> d == N+1 ? i : (d == N+2 ? j : Colon()), N+2)...]
        t2_comp = data(t2)[ntuple(d -> d == N+1 ? i : (d == N+2 ? j : Colon()), N+2)...]
        result_data .+= t1_comp .* t2_comp
    end
    
    ScalarField{T}(grid(t1), result_data)
end

"""
    outer_product(v1::VectorField, v2::VectorField) -> Tensor2Field
    v1 ⊗ v2

Compute outer product of two vector fields: (v1⊗v2)ᵢⱼ = v1ᵢv2ⱼ.
"""
function outer_product(v1::VectorField{T1,N,G}, v2::VectorField{T2,N,G}) where {T1,T2,N,G}
    compatible_fields(v1[1], v2[1]) || throw(ArgumentError("Incompatible grids"))
    
    T = promote_type(T1, T2)
    result_data = Array{T,N+2}(undef, size(grid(v1))..., N, N)
    
    for i in 1:N, j in 1:N
        v1_comp = data(v1)[ntuple(d -> d == N+1 ? i : Colon(), N+1)...]
        v2_comp = data(v2)[ntuple(d -> d == N+1 ? j : Colon(), N+1)...]
        result_data[ntuple(d -> d == N+1 ? i : (d == N+2 ? j : Colon()), N+2)...] = v1_comp .* v2_comp
    end
    
    Tensor2Field{T,N,G}(grid(v1), result_data)
end

# Tensor calculus operations

"""
    divergence(field::Tensor2Field, bc=nothing) -> VectorField

Compute tensor divergence: (∇⋅T)ᵢ = Σⱼ ∂Tᵢⱼ/∂xⱼ.
"""
function divergence(field::Tensor2Field{T,N}, bc=nothing) where {T,N}
    g = grid(field)
    result_data = Array{T,N+1}(undef, size(g)..., N)
    
    for i in 1:N
        component_data = zeros(T, size(g))
        for j in 1:N
            tensor_comp = data(field)[ntuple(d -> d == N+1 ? i : (d == N+2 ? j : Colon()), N+2)...]
            grad_comp = grid_gradient(g, tensor_comp, j, bc)
            component_data .+= grad_comp
        end
        result_data[ntuple(d -> d == N+1 ? i : Colon(), N+1)...] = component_data
    end
    
    VectorField{T,N,typeof(g)}(g, result_data)
end

# Mathematical operations

# Element-wise operations
function Base.:+(t1::Tensor2Field{T1,N,G}, t2::Tensor2Field{T2,N,G}) where {T1,T2,N,G}
    compatible_fields(t1[1,1], t2[1,1]) || throw(ArgumentError("Incompatible grids"))
    T = promote_type(T1, T2)
    Tensor2Field{T,N,G}(grid(t1), data(t1) .+ data(t2))
end

function Base.:-(t1::Tensor2Field{T1,N,G}, t2::Tensor2Field{T2,N,G}) where {T1,T2,N,G}
    compatible_fields(t1[1,1], t2[1,1]) || throw(ArgumentError("Incompatible grids"))
    T = promote_type(T1, T2)
    Tensor2Field{T,N,G}(grid(t1), data(t1) .- data(t2))
end

function Base.:*(field::Tensor2Field{T,N,G}, scalar::Number) where {T,N,G}
    S = promote_type(T, typeof(scalar))
    Tensor2Field{S,N,G}(grid(field), data(field) .* scalar)
end

Base.:*(scalar::Number, field::Tensor2Field) = field * scalar

function Base.:/(field::Tensor2Field{T,N,G}, scalar::Number) where {T,N,G}
    S = promote_type(T, typeof(scalar))
    Tensor2Field{S,N,G}(grid(field), data(field) ./ scalar)
end

function Base.:-(field::Tensor2Field{T,N,G}) where {T,N,G}
    Tensor2Field{T,N,G}(grid(field), -data(field))
end

# In-place operations
function Base.:+!(t1::Tensor2Field{T,N,G}, t2::Tensor2Field{T,N,G}) where {T,N,G}
    compatible_fields(t1[1,1], t2[1,1]) || throw(ArgumentError("Incompatible grids"))
    data(t1) .+= data(t2)
    t1
end

function Base.:-!(t1::Tensor2Field{T,N,G}, t2::Tensor2Field{T,N,G}) where {T,N,G}
    compatible_fields(t1[1,1], t2[1,1]) || throw(ArgumentError("Incompatible grids"))
    data(t1) .-= data(t2)
    t1
end

function Base.:*!(field::Tensor2Field, scalar::Number)
    data(field) .*= scalar
    field
end

function Base.:/!(field::Tensor2Field, scalar::Number)
    data(field) ./= scalar
    field
end

# Scalar field multiplication (component-wise)
function Base.:*(scalar_field::ScalarField{T1}, tensor_field::Tensor2Field{T2,N,G}) where {T1,T2,N,G}
    compatible_fields(scalar_field, tensor_field[1,1]) || throw(ArgumentError("Incompatible grids"))
    
    T = promote_type(T1, T2)
    result_data = Array{T,N+2}(undef, size(data(tensor_field)))
    
    scalar_data = data(scalar_field)
    for i in 1:N, j in 1:N
        result_data[ntuple(d -> d == N+1 ? i : (d == N+2 ? j : Colon()), N+2)...] = 
            scalar_data .* data(tensor_field)[ntuple(d -> d == N+1 ? i : (d == N+2 ? j : Colon()), N+2)...]
    end
    
    Tensor2Field{T,N,G}(grid(tensor_field), result_data)
end

Base.:*(tensor_field::Tensor2Field, scalar_field::ScalarField) = scalar_field * tensor_field

# Analysis functions

"""
    integrate(field::Tensor2Field) -> Matrix

Integrate each component of the tensor field over the domain.
"""
function integrate(field::Tensor2Field{T,N}) where {T,N}
    g = grid(field)
    results = Matrix{T}(undef, N, N)
    
    if has_volume_weighting(g)
        # For grids with special volume weighting (spherical, cylindrical)
        for i in 1:N, j in 1:N
            component_data = data(field)[ntuple(d -> d == N+1 ? i : (d == N+2 ? j : Colon()), N+2)...]
            results[i,j] = integrate(g, component_data)
        end
    else
        # Standard Cartesian integration
        vol_element = prod(spacing(g))
        for i in 1:N, j in 1:N
            component_data = data(field)[ntuple(d -> d == N+1 ? i : (d == N+2 ? j : Colon()), N+2)...]
            results[i,j] = sum(component_data) * vol_element
        end
    end
    
    results
end

"""
    frobenius_norm(field::Tensor2Field) -> ScalarField

Compute the Frobenius norm of the tensor field: √(Σᵢⱼ Tᵢⱼ²).
"""
function frobenius_norm(field::Tensor2Field{T,N}) where {T,N}
    result_data = zeros(T, size(grid(field)))
    
    for i in 1:N, j in 1:N
        component = data(field)[ntuple(d -> d == N+1 ? i : (d == N+2 ? j : Colon()), N+2)...]
        result_data .+= component .* component
    end
    
    result_data .= sqrt.(result_data)
    ScalarField{T}(grid(field), result_data)
end

# Alias for Frobenius norm
LinearAlgebra.norm(field::Tensor2Field) = frobenius_norm(field)

"""
    get_statistics(field::Tensor2Field) -> NamedTuple

Get statistical information about the tensor field components.
"""
function get_statistics(field::Tensor2Field{T,N}) where {T,N}
    component_stats = Matrix{NamedTuple}(undef, N, N)
    
    for i in 1:N, j in 1:N
        comp_data = data(field)[ntuple(d -> d == N+1 ? i : (d == N+2 ? j : Colon()), N+2)...]
        component_stats[i,j] = (
            min = minimum(comp_data),
            max = maximum(comp_data),
            mean = mean(comp_data),
            std = std(comp_data)
        )
    end
    
    frobenius_field = frobenius_norm(field)
    trace_field = trace(field)
    det_field = det(field)
    
    (
        components = component_stats,
        frobenius_norm = get_statistics(frobenius_field),
        trace = get_statistics(trace_field),
        determinant = get_statistics(det_field),
        integral = integrate(field)
    )
end

# Utility functions

"""
    to_scalar(field::Tensor2Field, method::Symbol=:frobenius_norm) -> ScalarField

Convert tensor field to scalar field using specified method.
Methods: :frobenius_norm, :trace, :determinant, :max_eigenvalue, :min_eigenvalue
"""
function to_scalar(field::Tensor2Field{T,N}, method::Symbol=:frobenius_norm) where {T,N}
    if method == :frobenius_norm
        return frobenius_norm(field)
    elseif method == :trace
        return trace(field)
    elseif method == :determinant
        return det(field)
    elseif method == :max_eigenvalue
        eigenvals = eigenvals(field)
        return eigenvals[1]  # First eigenvalue (sorted in descending order)
    elseif method == :min_eigenvalue
        eigenvals = eigenvals(field)
        return eigenvals[end]  # Last eigenvalue
    else
        throw(ArgumentError("Unknown method: $method"))
    end
end

# Display
function Base.show(io::IO, field::Tensor2Field{T,N}) where {T,N}
    print(io, "Tensor2Field{$T,$N}(")
    print(io, "grid=$(grid(field)), ")
    print(io, "size=$(size(field))")
    print(io, ")")
end

function Base.show(io::IO, ::MIME"text/plain", field::Tensor2Field{T,N}) where {T,N}
    println(io, "Tensor2Field{$T,$N}")
    println(io, "  Grid: $(grid(field))")
    println(io, "  Data size: $(size(data(field)))")
    println(io, "  Tensor rank: 2 ($N×$N)")
    
    # Show basic statistics
    stats = get_statistics(field)
    println(io, "  Frobenius norm: [$(stats.frobenius_norm.min), $(stats.frobenius_norm.max)]")
    println(io, "  Trace: [$(stats.trace.min), $(stats.trace.max)]")
    println(io, "  Determinant: [$(stats.determinant.min), $(stats.determinant.max)]")
end